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Abstract. 

We study coupled transport in the nonequilibrium stationary state of a model consisting 
of independent random walkers, moving along a one-dimensional channel, which carry a 
conserved energy-like quantity, with density and temperature gradients imposed by reservoirs 
at the ends of the channel. In our model, walkers interact with other walkers at the same site by 
sharing energy at each time step, but the amount of energy carried does not affect the motion 
of the walkers. We find that already in this simple model long-range correlations arise in the 
nonequilibrium stationary state which are similar to those observed in more realistic models 
of coupled transport. We derive an analytical expression for the source of these correlations, 
which we use to obtain semi-analytical results for the correlations themselves assuming a local- 
equilibrium hypothesis. These are in very good agreement with results from direct numerical 
simulations. 



PACS numbers: 66.10.cd, 05.60.Cd, 05.70.Ln, 05.40.Fb 

1. Introduction 

Two related outstanding problems at the heart of nonequilibrium statistical mechanics are 
the structure of the probability distribution function in the stationary state, and the derivation 
of macroscopic transport laws, such as Fourier's law of heat conduction, from microscopic 
dynamics, for systems which are maintained out of equilibrium by the imposition of 
thermodynamic fluxes [l] 12 . 

For certain classes of stochastic mass-transport equations, known as zero-range 
processes, in which the dynamics of mass leaving a site depends only on the occupation 
number at that site, the stationary-state distribution is known to factorise into the product of 
single-site distributions under certain conditions, which enables many analytical results to be 
obtained J4|. However, for more complicated models, this distribution no longer factorises. In 
this case, the appropriate characterisation of the stationary-state distribution becomes a central 
goal for the description of these systems. Furthermore, the fact that the distribution does not 
factorise implies the existence of spatial correlations between different sites, as has been 
discussed in many previous works |5j|6l. These spatial correlations in nonequilibrium states 
have been studied at a mescoscopic level using fluctuating hydrodynamics 0. Generically, 
they tend to be long-range, spanning the whole length of the system, rather than decaying 
exponentially as in equilibrium systems away from critical points. 
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From a microscopic point of view, the correlations arising in nonequilibrium stationary 
states have already been studied in many simple models, including a stochastic master 
equation describing heat flow [7], oscillators which exchange energy JSJ, lattice-gases with 
exclusion [9| 1 10], and lattice-gas cellular automata ifTTl . Also, exact results for all correlation 
functions were found using matrix technique for the symmetric simple exclusion process: see 
e.g. ifTZl for a recent review, where the relation of long-range correlations to a non-local 
free energy functional is also discussed. An approximation of the invariant measure using 
Gaussians in a suitably rotated coordinate system has recently been obtained lfl3l H4l : and 
related analytical results were previously found by Suarez et al. IfTTl . for the case where 
transport is by particles with exclusion, but with a single conserved quantity. 

In particular, long-range correlations in the so-called "random-halves" model lfl5l of 
coupled matter and heat transport were recently studied in 1161 . principally heuristically 
and numerically. The model we study in this paper can be considered to be a simplified 
version of the random-halves model, still containing two explicitly conserved quantities. The 
simplification enables us to obtain an explicit expression for the source of the correlations in 
the nonequilibrium stationary state of the system. 

Although our model suppresses much of the physical meaning of the second conserved 
quantity, in addition to mass, which in [ 16] can really be viewed as corresponding to energy, 
we emphasise that the structure of the spatial correlations that we observe for this energy-like 
quantity is remarkably similar to that found in ifTBl . 

In this paper, we study the equilibrium and nonequilibrium stationary states of the model. 
In regards to the nonequilibrium stationary state, we obtain the transport equations for the 
energy and mass, and we obtain the equation satisfied by the spatial energy correlations that 
arise in the model. This equation has a non-closed form. To close the hierarchy, we make 
a local equilibrium assumption, which enables an analytical evaluation of the form of the 
correlation source terms. We are then left with an approximate discrete Poisson equation with 
source terms for the correlations. We find very good agreement between the solution of this 
equation with numerical simulations of the system. 

2. Coupled transport model 

In this section we introduce the model of transport which we study. It is, perhaps, one of 
the simplest stochastic models exhibiting coupled transport. The transported quantities are 
particles (mass), and a second quantity, which is locally conserved, which the particles carry 
with them when they move. For simplicity of exposition, we refer to this second quantity as 
"energy", although we emphasise that it does not necesarily have the physical characteristics 
of an energy, since the motion of the walkers is independent of the value of the energy which 
they carry. 

Specifically, the model consists of independent random walkers moving on a one- 
dimensional chain of L sites. The system is open, and is in contact with particle and 
energy baths at each end of the chain, which at each time step supply or remove particles 
from the system with a given rate and energy distribution corresponding to their density and 
temperature, respectively. 

The walkers move synchronously in discrete time: at each time step, each walker 
independently attempts to jump to one of its two neighbouring sites, or remains at the same 
site. If a walker successfully jumps, then it carries with it an amount of energy s from the total 
amount of energy E at its previous site. 

After all particles have attempted their jumps, the total energy at each site, that is, the sum 
of the individual energies of the walkers at that site, is redistributed among all the particles 
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at that site, in a "random" (microcanonical) way, which we specify below. We thus have a 
complete (infinite) separation of time-scales: energy equilibration at each site is completed 
before the particles move again. This separation of time scales is, in part, what enables us 
to proceed with the analysis of the system. Further, it ensures that we can always use a 
local equilibrium hypothesis, in the sense that all thermodynamic quantities are always well 
defined at every site and that they are related to each other according to the usual (equilibrium) 
thermodynamic relations. 

As mentioned above, our model is related to the recently-introduced random-halves 
model lfT31 [161 . designed to model, rather faithfully, the Hamiltonian dynamics underlying 
the transport phenomena observed in |[T7l [T8l . In the random-halves model, each particle 
jumps to a neighbouring site with a rate which is proportional to the square root of its kinetic 
energy, and which is a factor 8 times the rate at which a particle exchanges a random fraction 
of its energy with a reservoir located at its current position. 

Taking the limit 8 — > in the random-halves model corresponds to our limit of infinite 
separation of time-scales, the other particles at that site acting as the reservoir. However, the 
random-halves model includes extra effects which our model cannot account for: for example, 
in the random-halves model, it is possible to have sites with particles at very low energy, and 
since the jump rate is energy dependent, these particles may remain a long time at that site 
unless other very energetic particles arrive there. Nonetheless, as we shall see, these kind of 
effects do not appear to affect the qualitative results on correlations. 

Although we do not consider it in this work, it should be noted that if we make the 
jump probabilities p and q very small (of order l/N, where N represents the number of 
particles in the system) we effectively recover single-particle motion, as in continuous-time 
dynamics, that is, on average, only one particle moves at each time step. Furthermore, in 
this case we could unambiguously consider jumping probabilities that are functions of the 
energy of the single moving particle, which could yield a model closer to that considered 
in fl6l . However, such modifications render the system intractable, and do not appear to 
be a necessary ingredient for the presence of long-range correlations in the nonequilibrium 
stationary states. 

2.1. Master equation 

We now proceed to specify the model precisely. We consider an arbitrary number of random 
walkers which can occupy sites on a finite one-dimensional chain of sites, labelled by 
i £ {1,...,L}. The system is open, and is in contact with particle and energy baths at sites 
and L+l. The baths have mean particle densities po and Pl+i, and are at temperatures 
To and Tl + \. This means that the number of particles available in each bath is drawn from a 
Poisson distribution with mean po and Pl+i, respectively, at each time step, and the energy E 
carried by each particle leaving a bath at temperature T has a Boltzmann distribution at that 
temperature, P(E) = je~ E l T . 

Let «, and £, be the number of particles and the total energy at site i at time f + 1, and 
m, and e, the corresponding quantities at time t. The walkers can jump to the right with 
probability p, jump to the left with probability q, or remain where they are with probability 
r := 1 — (p + q). The number of walkers which jump right from site i at a given time step is 
a random variable denoted if , and similarly, the number jumping left from that site is IJ . All 
walkers jump simultaneously. 

Each walker carries a certain amount of energy with it when it jumps. After each step, 
the new total energy Ej at a site i is distributed randomly among the n,- walkers at that site, 
according to a "microcanonical distribution". The total amount of energy carried by the 
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walkers which move from site i to the right is denoted by sf, and to the left by s t . 
The master equation describing the time evolution of this system is then given by 

P t+ \(n\,E\\n2,E2\ ■ ■ .;n L ,E L ) = 



£ £ / de, / dsf Pt(m l ,e l ;m 2 ,e 2 ;...;m L ,e L ) 

{<«.•} {/fr {e;} J{s ' } 

xl\8{n, - [ m + (/+ , + lf +1 ) - [if + If)]) (2 .1) 

i 

x [e, + (sU + sf +1 ) (st + sf)}) 

i 

x n p ( s i > s r | /, + , , «/, «*) x Il p ( z < + ' v I m '-) ■ 

i * 

The delta functions reflect the fact that the new occupation numbers and energies are obtained 
from the old ones by the movements at that time step. The conditional probabilities appearing 
in the last line of this equation denote the probability densities for the number of walkers and 
total energy moving left and right, and are given by 

P(/M» := (£)( ) pt (2.2) 

F(s+, sr\l+,ir,m hei ) := 



T(lt)T(li)T(mi-lf-ir 



trij — l 



The first is a trinomial distribution which gives the probability of moving exactly / ( + particles 
to the right and If to the left, out of the m, particles at site i. The second "multivariate 
beta distribution" is chosen to reflect the partitioning of the energy amongst the if, If and 
the remaining to, — / ( + — If particles, under the assumption that within each site the particles 
behave as a 2-dimensional ideal gas. 

The "ideality" of the gas at each site is manifested by the fact that the distribution can be 
written exactly in terms of products of appropriate phase space volumes, while the value of 
the exponents reflects the fact that the gas is 2-dimensional; details are given in the Appendix. 
This particular distribution was chosen because it yields slightly simpler expressions (than, 
say, 1 or 3 dimensional ideal gases) and is closer to the intrinsic 2-dimensional nature of 
various previous models for coupled transport. 

2.2. Equilibrium state 

It is known 1 19 1 that many non-interacting walkers, even when subjected to a density gradient, 
attain a stationary state which factorises: the probability of having occupation numbers 
n := (n\, . . . , n\f) is given by the following product of Poisson distributions at each site: 

L 

P(n)=l\P( ni ), (2.4) 



where 



e-P'O ni 

P(ni) = J 1 -, (2.5) 

nil 
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with p, the mean occupation number at site i. The p, satisfy a discrete diffusion equation, 
which in the stationary state becomes p, = pp,_i + qpi+i + rpi, and which can be solved in 
terms of the boundary conditions. 

Suppose now that there is no gradient of temperature imposed at the boundaries of our 
model for coupled transport, i.e. Tq = Ti + \ = T. Then it turns out the joint probability 
distribution of having energy Ei and n, particles at sites ; = 1,2... is also given by a factorised 
distribution: 

L 

P(E;n) =l\P(ni)P(Ei\ni), (2.6) 

;=1 

where the conditional probability of have energy £ at a site with n particles is given by 

f/npn-l -PE 

P(E\n) = ^-— , (2.7) 

i» 

where := l/T. (We take units such that the Boltzmann constant Icb = 1 throughout the 
paper.) This distribution can be interpreted as Q.(E,n)e~P E /Z(j3,n), where Q(E,n) is the 
volume of phase space accesible to a 2-dimensional ideal gas of n particles at total energy E, 
and Z(j3,n) is the partition function. 

The mean energy for this equilibrium distribution is (E\n) = n/fi, so that in equilibrium 
the mean energy at a site with mean concentration p is p//3 = pT. Since the distribution 
of energy is that of a system with temperature T, we can unambiguously identify with the 
inverse temperature. 

That the distribution factorises in equilibrium can be verified by assuming that the 
solution has a form as given in (12.61 l as an ansatz. It then transpires that the only way it 
can do so is if the temperature profile is flat. Hence, in the presence of a temperature gradient 
the joint distribution P(E;n) of all energies and positions does not factorise, and thus spatial 
correlations are present. 

3. Thermodynamics 

In this section we study the thermodynamic properties of the system. This is straightforward 
since, by construction, at each time step the system reaches a microcanonical equilibrium at 
each site i, characterized by the number of particles, «,, and the energy, , found at that site. 

Since we are assuming that at each site the particles are a 2-dimensional ideal gas, and 
accounting for the indistinguishability of the particles, the classical entropy at each site is 
given by 

S, =«,ln (^fj +nm, (3.1) 

where V is the volume (actually, the area) available for the gas at each site, which we take as 
unity (V = 1), and s, is a constant, in the sense that it is independent of the thermodynamic 
variables, though it may vary from one site to another (see Appendix). 

Having the fundamental relation (13.11 ). we can proceed to obtain the equations of state 
for the intensive variables in the entropy representation |20|: 

1 (dSi 



7; \dEjJ N . Ei 



and 

Ti ~ \dN, 



-lnff] + v, (3.3) 
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where v, := Si — 2 is a constant, which again may have different values at different sites. These 
expressions will be useful further on, in connection with the Onsager relations and the use of 
the local equilibrium hypthesis. 

3.1. Concentration and energy profiles 

We now consider the case in which the system is forced out of equilibrium by imposing 
concentration and/or temperature differences at the boundaries, that is, by imposing po ^ pi+i 
and/or Tq ^ Tl+\. If we do so, then the system will eventually attain a nonequilibrium 
stationary state, with well-defined concentration and mean energy profiles, p; and (£,), as 
a function of the position i in the system. Related profiles have been studied in detail for 
random-halves and other stochastic models in lfl"5ll2D . 

The transport equations can be easily derived by recalling that the total energy Ej at site 
i at time t + 1 is given by the energy e,- that was there at time t , plus the energy brought in by 
the walkers that arrived in that time step, minus the amount taken by the walkers that left: 

Ei = e x + (*+ , + sr +l ) - (s+ + sr ) . (3.4) 

Taking means, we obtain 

(Ei) = (e i )+p(e i _i)+q{e i+l )-p{e i )-q{e i ) = + q(e i+i ) + r{e t ) . (3.5) 

A similar equation holds for particle transport. In the stationary state, (Ei) = (ei), and hence 
the stationary profiles satisfy 

(tii) = p (ni-x) + r (m) + q (n i+ i) ; (3.6) 
(E t ) =p(E i _ l )+r(E i )+q(E i+l ). (3.7) 

We denote by p, := («,•) the stationary mean occupation number at site i, and by 
Ti := (Ei) I pi the local temperature there. 

3.2. Thermodynamic fluxes and forces 

The mean energy and particle fluxes between sites i and i + 1 are given by 

/„ = p (Ei) - q (E i+ i) = pPiTj - qp i+l T i+u (3 . 8a) 

Jp= ppi-qpi+i. (3.8£) 

To obtain the continuum (diffusive) limit, we first express p,, 7], p i+ i and 7} + i as Taylor series 
around position x = (/+ 1/2) 8x, where 8x is the distance between neighbouring sites on 
the chain. Next, we transform p — > c8x, J u —> j u 8t and J p —> j p 8t, where 8t is the time 
interval between succesive steps, so the quantities c, j u and j p are a proper density and fluxes, 
respectively. These operations yield 

j u St = {p-q)cT8x-\{p + q){8x) 2 V{cT) + ff{8x 3 ) (3.9a) 

jpSt = (p-q)cSx -\{p + q){8x) 2 Vc + 0{8x^). (3.9b) 

The continuum limit is achieved by dividing through by St and taking the limit in which 
8t, 8x and p — q tend to 0, in such a way that the ratios (8x) 2 /8t and (p — q)/ Sx remain finite. 
Thus, we can define the drift velocity 

v:= lim (p-q)Sx/St (3.10) 



and the diffusion constant 

2 



D :=l lim (p + q) (8x) 2 /8t, (3.11) 
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in terms of which the above equations become 

j„ = -DV(cT) + vcT (3.12a) 
jp = -DVc + vc. (3.12b) 

These should be compared with 

j u =L 11 V(l/7')+L 12 V(- A t/r) (3.13a) 
j p =L 2l V(l/T)+L 22 V(-fi/T) (3.13b) 

from the theory of linear thermodynamics [20, 22]. Using d3.2| ) and j3.3K we obtain 

j u =-L n (Vr)/r 2 +L 12 [(Vr)/r-(V c )/c + Vv] (3.14a) 
jp = -La {VT)/T 2 +L 22 [(VT)/r - (Vc)/c + Vv] . (3.14&) 

Setting Vr = Vc = 0, we find 

L l2 Vv = vcT and L22VV = vc. (3.15) 

If instead we set VT = Vv = 0, then 

L l2 /c = DT and L 22 /c = D. (3.16) 

Finally, if we set Vc = Vv = 0, then 

L u /T 2 -L 12 /T=Dc and L 22 /T -L 21 /T 2 = 0. (3.17) 

From these equations, we obtain: 

L u =2DT 2 c; L l2 = L 2l = DTc; (3.18) 

La=Dc; Vv = v/D. (3.19) 

Thus the Onsager reciprocal relations Ell are satisfied, and v, is determined as an external 
potential due to the overall current generated by the bias. 

While these results are satisfactory, it should be noted that we made a rather cavalier use 
of (13.2b and (13.31 >. namely, we identified the temperature at site i as the quantity (£,•) / («,}, 
whereas ( 13.2b tells us that the local temperature is actually the stochastic variable 7j- = Ei/nf, 
furthermore, we substituted the remaining in i3.3\ by c, 8x — p, = («,•). These substitutions 
are, of course, not generally valid; however, they can be justified when the fluctuations in 
energy and number of particles are small compared to their mean values. 



4. Spatial correlations of the energy 

We now turn to the main consideration of the paper, the origin of the spatial correlations 
between the values of energy at different sites, which develop due to the imposition of a 
temperature gradient. To this end, we denote by Qj := {EiEj) — (£,■) (Ej) the stationary- 
state energy correlations between sites i and j. To simplify the notation, we use the difference 
operators A™ and A/,-) which act on functions of two variables Qj as 

[ A (o c ] 'j '■ = pd- 1 j + rC 'j + <i c i+ hj< ( 4 - 1 ) 

[A(j)C];j := pQj-i + rQj + qQj+i . (4.2) 
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4.1. Exact equation for stationary-state energy correlations 

Using the above notation, it follows from the previous section that the evolution equation for 
the average energy is (Ej) — Am (e,), and thus 

<fr)(Ej)=\ i) A U) [(e i )(ej)], (4.3) 

so that this part of the correlations Qi factorises. 

It remains to evaluate (EiEA. To do so, we rewrite it quantity in terms of the energies e,- 
and ej at sites i and j before the move, and the amounts of energy moving in each direction 
from each site: 

(EiEj) = ([ei+(4-i+'Hi)-(4+'T)] ■ h + (4-i+' s 7+i)-(4 +i 7)])- < 4 - 4 > 

We expand the product and consider the resulting terms, which are means of products of 
two random variables, of the form (^4-i4+i) ' According to the master equation d2.lt . 
these random variables are independent if their indices are different, giving, for example, 
(^l^j+i) = (4i) if ' — 1 7^ ./' + 1- in particular, this is the case for every pair of 

products provided |/ — j\ > 2. 

If, on the other hand, \i — j\ < 2, then there are terms in the product for which the indices 
are the same: for example, if j = i+ 1, then s~ i+l = sj. In this case, the mean of the product is 
no longer the product of the means, and we must calculate it explicitly. For example, we have 

l + n l + (l + + l)e 2 

{st\eult,m) = {4 2 \*i4,m) = 1 [ ) V • (4.5) 

x ' nit x ii nii{mi + 1) 

We must then average the expressions over the trinomial distribution for the if and If given 
m. Note that the right-hand side of the second equation is not the square of the first equation - 
a correction term has arisen. These corrections are what eventually give rise to the long-range 
energy correlations. 

We finally obtain, after some messy algebra, which we confirmed via a computer algebra 
package, the following exact equation for the spatial correlations Cu in the stationary state: 

C i , / = A (I) A (; - ) C, J +2A (7 , (4.6) 

with Ay a symmetric matrix given by 

>(1 -/?)«■,_!+ r(l - r) K t + g(l - q) K i+1 if j = i; 

, -prK; — qrK: 4. 1 if j = i+ 1; 

Xa=i ' H 1+1 J (4.7) 

1 -pqKi+i rfj = i + 2; 

otherwise, 
where we have defined 



K, := ( —j— ) . (4.8) 
\ nij + 1 / 

Equation (14.61 > is essentially a discrete Poisson equation, with source terms 2A;;. 

The boundary conditions are Qj = whenever i or j is equal to or L+ 1, since the 
stochastic reservoirs at positions and L + 1 are independent of all other quantities in the 
system (and of each other). The exceptions to this are Cq.q and Cl+i.l+\, which are given by 
the variances of the distributions in the reservoirs. 

The above equations may be simplified by introducing 

*y:=Cij-2fyjQ, (4.9) 
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where 5, 7 is the Kronecker delta, which is 1 when i = j and otherwise. Substituting this 
expression in ( 14.61 ) gives that gy satisfies the following simpler equation: 

g u A, A l)giJ • 2u : S, r (4.10) 

with source terms 

fii := pKj_! + (r-l)Kj. (4.11) 

Note that the only source terms in ( 14.101 ) are now on the diagonal. The boundary conditions 
are g, j = when i = 0, i = L + 1, j = or j = L + 1 . We thus have a discrete Poisson equation 
in a square, with zero boundary conditions and a line source term on the diagonal. 

We can test this equation in the simplest case: that in which there is no energy 
(temperature) gradient. In this case, the reservoirs are at the same temperature j3 := j3o = Pl+i> 
so that in fact the temperature is constant throughout the system, p\- = j3 for all i. Under these 
conditions, we know that the energy distribution factorises, hence there must be no energy 
cross-correlations. Indeed, in this case we can evaluate JQ exactly to obtain fq = p,/j3 2 , and 
( 14.111 ) then gives 

Hi = iPPi-i +qpi+\ + (r- I) pi\ = 0, (4.12) 

since the p, satisfy precisely this discrete equation. Thus, in the absence of a temperature 
gradient, the gjj satisfy g(j = A/{\Ar:\gij for all i and j, with no source terms. The zero 
boundary conditions then imply that gu is identically zero. 

Substituting this result back into ( 14.9b , we obtain in this constant temperature case 

Qj =25,^=2^%. (4.13) 

The term 25, ; fq can thus be regarded as the contribution to the energy correlation matrix 
which arises simply because C,-., necessarily has a non-zero on-site value, given by 

Q, = (E? )-(Ei) 2 =2 K , = ^. (4. 14) 

Referring back to the definition (14.91 ) of we see that this quantity can thus be 
viewed as containing the long-range part of the correlations, resulting from the imposition 
of temperature gradients. This is similar to results of previous work in the case of a single 
transported quantity J7] QT) . 

We remark that the physical meaning of the terms k - ,- = (ef/(m.i +1)), which form the 
source terms of the long-range correlations, and thus in some sense are what gives rise to these 
correlations, is not very clear. 



4.2. Local thermodynamic equilibrium approximation 

The previous calculation is exact; however, to make further progress, we must make an 
approximation in order to evaluate the terms jq appearing in the expression for the source 
Hi of the long-range part of the correlations when the system is in a nonequilibrium stationary 
state. To do so, we will assume that local thermodynamic equilibrium is attained at each 
site. By this we mean the assumption that the marginal distribution of the energy at each site 
i is given by P(Ei\ni), with the distribution (12.71 ) which is found at equilibrium. This is an 
uncontrolled approximation; however, we will see later that it is in very good agreement with 
the numerical results. Note that jc, involves only data at site i, and thus indeed depends only 
on the marginal distribution at that site. Such a local thermodynamic equilibrium assumption 
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z -piTf. (4.15) 



has recently been proved correct for the random-halves model, in the limit when the number 
of sites goes to °°, so that the temperature gradient goes to zero |23l . 

Under the hypothesis of local thermodynamic equilibrium, we can use (12. 7t to calculate 

Kj = (Ef/{mi + 1)), obtaining 

We will use this approximate analytical form for (Q in the remainder of the analytical 
development. 

From the above discussion, we see that the contributions to the onsite correlations Cjj 
split into two parts: 

C lti =g iti + 2rq. (4.16) 

We can regard 2fC, as a pure local contribution, and g, - j as the onsite part of the long-range 
contribution. Within the local equilibrium approximation, we then obtain 

g u = Cu-2 Pi T?. (4.17) 

The equations for the profiles of mean concentration and mean energy can be written as 
follows: 

ppi-r+qpi+y + ir-^pt =0; (4.18) 

PPi-iTi-! +qp i+1 T M + (r-l)PiTi = 0. (4.19) 

For brevity, we introduce the linear operator Jz^ [f] : = pfi- \ + qfj+ \ + (r — 1)/; = (Ain — 1 )/. 
The equations for the profiles then become ££i\p\ =0 and Jzfi[pT] = 0. 

The correlation source is /i, = Jzf/[fc] = J^j{pT 2 ], where the latter equality again assumes 
the local thermodynamic equilibrium approximation. Substituting J^ [p] = and J^ [pr] = 
into the expression for ^[pT 2 ], we obtain that the source term in the local thermodynamic 
equilibrium approximation is given by 

Mi = ^i[pT 2 } = \{T i+x - 7}_i) \ppi-i{Ti - Ti-x)+qp i+ i{T i+l - 7]-)] . (4.20) 

In the continuum diffusion limit, we have 

&\f\~[Df- V f\8t, (4.21) 

with v, D are defined in equations ( 13.1 0b and ( 13.111 ), respectively, and p = c8x as before. 
Then the source of correlations fi (x) reduces to 

H(x) ~ ^[pT 2 ] ~ 2Dc(x)[T'(x)] 2 StSx, (4.22) 

as can also be verified directly from the continuum expressions. We remark that precisely 
a quadratic dependence on the local temperature gradient of short-range energy correlations 
was found numerically for the random-halves model 1 16 1. 

It should be noted that only the energy has long-range correlations: a calculation 
similar to the above shows that the density correlations («;«;) and the density-energy cross- 
correlations (Eitij) are both diagonal. 



5. Numerical results 



In this section, we present comparisons of the energy correlations as obtained from direct 
simulations of the microscopic random-walk dynamics, with the approximate analytical 
results derived in the previous section. 
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Figure 1. (a) Long-range part of the correlations, gf< = Qi — 2&/1Q, for an imposed 
temperature gradient and flat density profile, with L = 41, p = 10, To = 50, 7/, + i = 10. 
Shown are the direct numerical results, using the numerical values of Cu and IC,-, and the 
semi-analytical approximation obtained by solving numerically the discrete Poisson equation 
44. 101 using the local equilibrium assumption K, ~ p,//3 2 . In this and the following figures, 
each separate curve shows the stationary-state energy correlations g, ■ j for a given position i as a 
function of j. (b) Rescaled correlations Lg t j in the absence of a density gradient, as a function 
of j/L, i = (L + 1 )/2 and i = (L + 3) /4, for the same parameter values as in (a) and different 
system sizes L. Onsite correlations g l j are not shown. The exact result in the continuum limit 
is shown as a solid line. 



5.7. Numerical method 

The boundary conditions in the numerical simulations are as follows. At each time step, the 
number of particles «o at the left bath is chosen from a Poisson distribution with mean po, 
and each of those particles is assigned an energy E with probability l/Toe~ E / T ° . The same 
is done at the right bath with appropriate temperature and density. At each site, the energy 
of the particles is assigned via the microcanonical redistribution mechanism, and then each 
particle jumps to a neighbouring site with the correct probabilities. Means and correlations 
are determined by time averaging, after discarding a preliminary equilibration period. 

The correlations from the direct numerical simulations are compared to "semi-analytical" 
results obtained by solving the discrete Poisson equations (14.91 ) for the long-range part g of 
the correlations, using the local equilibrium approximation ( I4.15l l for the terms Kj. A similar 
numerical solution of the algebraic equations was recently employed in |[T3l . 

5.2. Temperature gradient in absence of concentration gradient 

The simplest case with non-trivial correlations is to impose a linear temperature gradient 
between two baths with temperatures Tq ^ TL+i, but with a flat concentration gradient, i.e. 
Po = Pl+\ = P for all U and without bias in the dynamics (p ~ q). In this case, the density 
profile is flat throughout the system, p, = p for all i. The profile of mean energy is linear: 

(E,) = (E Q ) + -±- x ((E L+l ) - (E )) • (5.1) 
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Figure 2. Difference AKj := (f, — p 7} 2 of the value of IC, obtained from the direct numerical 
simulation and the analytical expression obtained with the local equilibrium assumption, for 
two values of system size L. The differences A/i, between the corresponding results for /J, 
obtained by differenciation are also shown. 



Identifying as usual 7J := (£,) /pi, we can conclude that there is also a linear temperature 
profile under thse conditions; this is correctly obtained in simulations (not shown). 

Figure [T(a)| shows the long-range part of the energy correlations, gij, in this case. 
Numerical results are compared to a numerical solution of the algebraic discrete Poisson 
equation ( 14.9b . In order to carry out this numerical solution, the source terms ff, were assumed 
to take their local equilibrium value p as described above. Despite this, we find very good 
agreement between the numerical results obtained from direct simulation and the numerical 
solution of the discrete diffusion equation. This holds everywhere, including for the onsite 
contribution of g,- ,-. 

Nonetheless, the agreement between the numerical and semi-analytical results is affected 
by the fact that the local thermodynamic equilibrium approximation is not strictly correct. 
As discussed in the introduction, the structure of the out-of-equilibrium measure is an open 
problem. However, here we can obtain an indication of the error in the local thermodynamic 
equilibrium approximation by comparing the value of kj = (e?/(m, + 1)) obtained from a 
direct numerical average to the analytical value p,/ /3? obtained from the local thermodynamic 
equilibrium assumption. This difference is shown in figure|2]for two different values of system 
size L. We see that the marginal distribution is not quite given by the local equilibrium 
approximation, but that it gets closer as L increases, in agreement with the rigorous results 

of JH. 

For the structure of the correlations gu away from the diagonal terms where i = j, 
the important quantities are the sources which are given by differences of the jq as in 
( 14.1 II ). The difference between the /I, calculated by differentiating the numerically-obtained 
fQ, and those obtained by differentiating the local equilibrium expression, are also shown in 
figure |2 They are very close to 0, which is the reason for the excellent agreement between 
the numerical and semi-analytical results for the correlations. 

In fact, this case (absence of concentration gradient) is simple enough to solve explicitly 
in the continuum limit. Taking g(i. j) — > G(x,y) 8x 8y, equation ( 14.1 Oi l can be rewritten in the 




Figure 3. Long-range part of the correlations, g, with imposed temperature and density 
gradients. The parameters are L = 41, Tq = 50, and Tj, + i = 10. The bath densities are po = 10 
and Pl+i = 20 in (a), and are reversed in (b). The direct numerical and semi-analytical results 
again agree very well. 



continuum limit as: 

dG(x,y) dG(x,y) 



:4c[T'(x)] 2 S(x-y), (5.2) 



dx 2 dy 2 

with boundary conditions G(x = 0,y) — G(x = L,y) = G(x,y = 0) = G(x,y = L) = 0. The 
solution of this equation is readily found to be 

G(x,y) = ^j^-l . (5.3) 

[ y(L — x), ify<x 

This result is similar to those of refs. 0QT), but with the difference that the concentration c 
now appears explicitly in the result. 

As was pointed out in fl6l . the correlations for a system of size L decay as l/L if 
the boundary conditions are fixed (i.e. the values of the density and temperatures at the 
boundary are the same for diferent values of the system size L). Figure [T(b")| shows the rescaled 
correlations Lg, ;j for different system sizes compared to d5.31 l. In the figure we have thus 
scaled space to the interval [0, 1] and rescaled the correlations by multiplying them by L. The 
various curves indeed converge to the limiting (continuum) form as L — > °°. Note that the 
apparent l/L scaling arises from the 8x term in the passage to the continuum limit: fixing the 
total system size and doubling the number of sites corresponds to halving 8x. 



5.3. Combined temperature gradient and concentration gradient 

We now consider the effect of imposing both energy and concentration gradients, although 
still without bias in the motion (p = q). The profiles of p,- and of (£;) = p,7] are now both 
linear, so that 7} = (£,-) / p,- is a ratio of two linear functions, but is no longer itself linear. 

Figure[3]shows the numerical and semi-analytical results in this case. We find a skewing 
effect on the correlations, which is again in excellent agreement with the numerical results. 




0.2 0.4 0.6 0.8 1 



j/L 

Figure 4. Scaling of correlations Lgu, obtained from the semi-analytical solution, for different 
system sizes L, with parameters as in figure[5Ja), for i = (L + 1 ) /2 and i = (L + 3 ) /4. 



The results are qualitatively very similar to those in fl6l . despite the differences in the nature 
of the models discussed in the introduction. Figure [4] shows the scaling of the correlations 
(obtained from the semi-analytical results) with system size. Again they converge to a 
continuum limit, corresponding to the solution of the continuous diffusion equation with 
sources given by ( 14.221) . 

5.4. Effect of bias (p ^ q) 

Upon introducing a bias in the directionality of the walkers' jumps, that is by putting p ^ q, 
we obtain mean density and energy profiles which are no longer linear. Rather, they are given 

by m 

1 - a' 

1 — O!^ 1 

1 - a' 

(Ei) = (Eo) + l aL+l ({E L+l ) - (Eb» , (5.5) 

where a := p/q and the quantities Pl+i, Po, (El+i) and (Eq) are fixed by the boundary 
conditions. 

Figure [5(a)| shows the comparison between the numerically-obtained correlations and the 
semi-analytical solution, for a situation with a fiat density profile (i.e. Pl+\ — PoX an imposed 
temperature gradient (which for constant density means that (£l+i} (Eq)) and a bias to the 
left. Again we find excellent agreement between the numerical results and the semi-analytical 
results. Note, however, that for values of i around 30 and larger, the correlation function gij is 
very small for all values of j (j ^ ;')■ Hence, the energies at sites corresponding to large enough 
values of i and j are essentially uncorrelated. The reason for this is that, in the presence of 
a bias, the source of the correlations given in eq. (14.22b decays exponentially as i increases. 
Thus, in this case, since the bias is to the left, we expect the source of the correlations to 
be appreciable only up to distances of a few times the decay length A ~ 1 / 1 ln(oc) | (see eq. 
( 15.5b ) from the left boundary, and the effect of the right boundary to become negligible if the 
size of the system L >> X. Indeed, figure [5(b)l compares the correlations for two different 
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Figure 5. (a) Comparison of numerical and semi-analytical correlations g, <j for system size 
L = 41, with parameters p = 10, To = 50, 1 = 5, as a function of j, for ;' = 3, 4, 5, 8, 11, 
16, 21 and 30 from top to bottom. There is a bias in the dynamics, with p = 0.35 and q = 0.4. 
(b) Comparison of the correlations for L = 41 and L = 81, with the same parameters as in (a). 
The numerical value of the correlations are the same for not too large. The lines are shown 
as a guide for the eye 
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Figure 6. Rescaled semi-analytical correlations Lgjj for biased dynamics as a function of 
j/L, for system sizes L = 41,L = 81 and L = 161 and different i. In order to have a well- 
defined continuum limit, q — p is halved when the system size is doubled, with q = 0.4 fixed. 
Parameters are p = 10, To = 50, 7i + i = 5. 



system sizes with the same parameters. Clearly, the numerical values of the correlations are 
essentially independent of system size, and are negligible over a large part of the system. 
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Nonetheless, as was shown in section 13.21 it is possible to obtain a well-defined 
continuum limit for the energy correlations in the case of biased dynamics, provided the 
amount of bias changes in the correct way: the bias q — p must be halved when the system 
size doubles, in addition to performing the same linear rescaling of space as in the other cases. 
Figure [6] confirms convergence to the continuum limit under these conditions. 

6. Conclusions 

By studying an extremely simple model of coupled transport of mass and a second conserved 
quantity, which we called energy, we have shown that the ubiquitous long-range correlations 
in this energy, whose transport depends on the motion of the mass, are present in the 
nonequilibrium stationary state, even though this quantity is completely passive. We were able 
to write down the equation describing the long-range spatial correlations for the energy in the 
system, and found that the structure of these correlations is remarkably similar to those found 
in more realistic models, as well as to the results of studies using fluctuating hydrodynamics. 
We thus conclude that the origin of these long-range correlations is already present in this 
simple model, and that a study of such models can go at least part of the way to explaining and 
quantifying the origin and structure of correlations in nonequilibrium systems. Nevertheless, 
we hope to be able to extend the methods and results to cases where the particle motion is 
modified by the energy carried by the particles. 
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Appendix A. Derivation of energy partitioning distribution and entropy at a site 

The distribution for the microcanonical pardoning of energy among the particles at each site 
can be calculated as follows. First, given that the particles at each site are assumed to be in 
a state of microcanonical equilibrium, we can determine the probability that i ( - and IJ out of 
the rrii particles at site i have combined energies s~f and sj respectively, when the total energy 
at the site is e,-, as the quotient of the number of states consistent with these requirements 
over the total number of states available to the system. These numbers are proportional to the 
corresponding structure functions [24], given by 



where i/^(p,q) is the A^-particle Hamiltonian describing the dynamics at each site, and the 
integration is carried out over the phase space of the N particles. If we assume that the particles 
are an ideal gas, then 




(A.l) 
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where m is the mass of the particles, and the required probability is given by 

W{st,sr\l+,ir, muei )= (A.3) 
co (sf,lf, V ) co (sj , If , V ) co (gj - - jr , m , - /+ - it , y ) 
©(e,-,mi,y) 

For a (i-dimensional ideal gas, we have 

^• V » = U « < (A ' 5) 

where V is the volume accesible to the particles at each site, which we take as unity, and T(-) 
is the gamma function. Using this expression in jA.3i yields 

F(sf,s-r\l+,ir, mi , ei )= (A.6) 

T{ mi d/2) 



T{l+d/2)T{i;d/2)T{[mi - if - ir]d/2) 

l P . ,+ „-\[m,-l?-ir]d/2-l 

( +y+rf/ 2 -i f -yrrf/2-i C e » s i ~ s i > [ ' ' 



mid 12— \ 



(A.7) 
(A.8) 



which, when d = 2, simplifies somewhat, giving 

P(*+ | /+,/,-, m/,«i) = (A.9) 

nm) fj+ ^-i (j -^-i CgiZitz£gll^ (A 10) 

r(/+)r(/r ) r(m i -/ ( + -/r)^ ) ■ (A ' 10) 

Similarly, the classical entropy of the 2D gas at each site is given by |20| 

S(e,-,m,-,y) = ln[©(e,-,mi,7)/w,!] +m ! a i - m,-ln[Ve,-/m 2 ] +m,-j,-, (A.11) 

taking the Boltzmann constant = 1, where s, and a, are constants (i.e. they are independent 
of e,, m, and V) which can vary from site to site. 
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